home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / QFLOAT / QPOW.C < prev    next >
C/C++ Source or Header  |  1996-02-25  |  2KB  |  155 lines
  1. /*                        qpow.c    */
  2. /*  power function: z = x**y */
  3.  
  4. #include "qhead.h"
  5. #include "mconf.h"
  6.  
  7. extern QELT qone[];
  8. int qpowi(), mtherr();
  9.  
  10. int qpow( x, y, z )
  11. QELT *x, *y, *z;
  12. {
  13. QELT w[NQ];
  14. long li;
  15.  
  16. qfloor( y, w );
  17. if( qcmp(y,w) == 0 )
  18.     {
  19.     qifrac( y, &li, w );
  20.     if( li < 0 )
  21.         li = -li;
  22.      if( li < 16384L )        /* 02-20-96 per Steve Moshier */
  23.         {
  24.         qpowi( x, y, z );
  25.         return 0;
  26.         }
  27.     }
  28. /* z = exp( y * log(x) ) */
  29.  
  30. qlog( x, w );
  31. qmul( y, w, w );
  32. qexp( w, z );
  33. return 0;
  34. }
  35.  
  36.  
  37. /* y is integer valued. */
  38.  
  39. int qpowi( x, y, z )
  40. QELT x[], y[], z[];
  41. {
  42. QELT w[NQ];
  43. long li, lx;
  44. int signx, signy;
  45.  
  46. qifrac( y, &li, w );
  47. if( li < 0 )
  48.     lx = -li;
  49. else
  50.     lx = li;
  51.  
  52. if( lx == 0x7fffffff )
  53.     {
  54.     qpow( x, y, z );
  55.     return 0;
  56.     }
  57.  
  58.  
  59. if( x[1] == 0 )
  60.     {
  61.     if( li == 0 )
  62.         {
  63.         qmov( qone, z );
  64.         return 0;
  65.         }
  66.     else if( li < 0 )
  67.         {
  68.         qinfin( z );
  69.         return 0;
  70.         }
  71.     else
  72.         {
  73.         qclear( z );
  74.         return 0;
  75.         }
  76.     }
  77.  
  78. if( li == 0L )
  79.     {
  80.     qmov( qone, z );
  81.     return 0;
  82.     }
  83.  
  84. qmov( x, w );
  85. signx = w[0];
  86. w[0] = 0;
  87.  
  88.  
  89. if( li < 0 )
  90.     {
  91.     li = -li;
  92.     signy = -1;
  93.     }
  94. else
  95.     signy = 0;
  96.  
  97.  
  98. /* First bit of the power */
  99. if( li & 1 )
  100.     {
  101.     qmov( w, z );
  102.     }
  103. else
  104.     {
  105.     qmov( qone, z );
  106.     signx = 0;
  107.     }
  108.  
  109. /* Overflow detection */
  110. /*
  111. lx = signy * li * (long)w[1];
  112. if( lx > (long)MAXEXP )
  113.     {
  114.     qinfin( z );
  115.     mtherr( "qpowi", OVERFLOW );
  116.     goto done;
  117.     }
  118. if( lx <= 0 )
  119.     {
  120.     qclear( z );
  121.     return 0;
  122.     }
  123. */
  124.  
  125. li >>= 1;
  126. while( li != 0L )
  127.     {
  128.     qmul( w, w, w );    /* arg to the 2-to-the-kth power */
  129.     if( li & 1L )    /* if that bit is set, then include in product */
  130.         qmul( w, z, z );
  131.     li >>= 1;
  132.     }
  133.  
  134. /*
  135. done:
  136. */
  137.  
  138. if( signx )
  139.     qneg( z ); /* odd power of negative number */
  140.  
  141. if( signy )
  142.     {
  143.     if( z[1] != 0 )
  144.         {
  145.         qdiv( z, qone, z );
  146.         }
  147.     else
  148.         {
  149.         qinfin( z );
  150.         mtherr( "qpowi", OVERFLOW );
  151.         }
  152.     }
  153. return 0;
  154. }
  155.